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1 Introduction 

There is a tendency to write the equations of general relativity as a first order 
symmetric system of time dependent partial differential equations. However, for 
. numerical reasons, it might be advantageous to use a second order formulation 

^ \ like one obtained from the ADM equations. Unfortunately, the type of the ADM 

equations is not well understood and therefore we shall discuss, in the next 
section, the concept of wellposedness. We have to distinguish between weakly 
. and strongly hyperbolic systems. Strongly hyperbolic systems are well behaved 

\ even if we add lower order terms. In contrast; for every weakly hyperbolic system 

we can find lower order terms which make the problem totally illposed. Thus, 
for weakly hyperbolic systems, there is only a restricted class of lower order 
O . perturbations which do not destroy the wellposedness. To identify that class can 

be very difficult, especially for nonlinear perturbations. In Section 3 we will show 
5-H . that the ADM equations, linearized around flat with constant lapse function and 

shift vector, are only weakly hyperbolic. However, we can use the trace of the 
metric as a lapse function to make the equations into a strongly second order 
k> " hyperbolic system. 

. Using simple models we shall, in Section 4, demonstrate that approximations 

' of second order equations have better accuracy properties than the corresponding 

approximations of first order equations. Also, we avoid spurious waves which 
travel against the characteristic direction. 

In the last section we discuss some difficulties connected with the preservation 
of constraints. 

2 Well Posed Problems 
2.1 First Order Systems 

Consider the Cauchy problem for a first order system with constant coefficients 

= E ^^^^^ =■■ ^(^)"' = 7^.^ (1) 
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u{t = 0) = f{x), X = (xi, . . . , Xs), — OO < Xj < CO. 

We construct simple wave solutions 

u{x,t) ^ e''''^'''^^u{uj,t), Lj = {uji, . . . ,u!s) real 

and obtain 

s 

^ i\uj\^AjUjjU =: i\uj\P{uj')u, (2) 

Definition 1. We call the problem (|l|) strongly hyperbolic if for every uj' the 
eigenvalues of the symbol P{oj') are real and there is a complete set of uniformly 
(in uj') linearly independent eigenvectors. 

Examples of strongly hyperbolic systems arc symmetric systems where Aj = 
The solutions of strongly hyperbolic systems satisfy an energy estimate 

||«(-,t)||'<^V"*||tt(.,0)ir. (3) 

Here K, a are universal constants which do not depend on the initial data 
u(a;,0) = f{x). The norms are L2 norms. For systems with constant co- 
efficients a = 0. 

Strong hyperbolicity and the existence of an energy estimate are equivalent, 
we hav^ 

Theorem 1. The solutions of (|l|) satisfy an energy estimate of the type (H) if 
and only if the system is strongly hyperbolic. 

The most important property of strongly hyperbolic systems is that we can 
add lower order terms and an estimate of type (||) is still valid. We have 

Theorem 2. Let ([^) be strongly hyperbolic. Then the solutions of 

wt = P{D)w + Bw (4) 
satisfy an estimate of type (^) . Here B is any bounded operator. 

Definition 2. We call the problem (|l]) weakly hyperbolic if the eigenvalues of 
P{u!') are real. 

In this definition we do not require that there is a complete set of eigenvectors. 
An example in dimension one is given by 

= ( Q M Mi; =: Au^ 



^ First order theory is well known, we refer to 
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Simple wave solutions for this system have the form 

■u(a;,t) =e''^^V'^^S(w,0) = + ioj (j^^^^ e''^(^+*)u(a;, 0) 

Thus there is no exponential growth but there is the polynomial growth in ivt. 
This is typically the case for weakly hyperbolic systems. One can prove 

Theorem 3. For weakly hyperbolic systern,s the growth of simple wave solutions 
is at most of the order 0(1 + Iwil""^), where n is the number of components of 
u. 

The real difficulty with weakly hyperbolic systems is that lower order terms 
can make them exponentially ill posed. For example, consider 



"*=ioij"^+lio'"- 

Making the simple wave ansatz 

u{x,t)=e''^='u{uj,t) 

we obtain 

^ f iu! iuj\ ^ . ^ 
Ut=[.. \u=: Au. 

The eigenvalues A of ^ are given by 

/r- \/2 / 

A = ia;±viw, i.e. ReA = ±^-i/|a;|. 

Therefore the perturbed problem is exponentially ill posed. 

The same result holds if we consider the variable coefficient problem 

There are no problems to generalize the results to variable coefficients and 
quasi-linear systems if the system is pointwise strongly hyperbolic. 

Unfortunately, in applications one can be confronted with systems which are 
weakly hyperbolic. In this case one has to carefully study the influence of lower 
order terms. For example, trivially, 

ut + ('W^)x + {v'^)x = -au, vt + {v^)x = -av 

is well behaved (a > sufficiently large so that no shocks form). We can solve 
the second equation to obtain v which becomes a governing function in the first 
equation. 
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2.2 Second Order Systems 

We consider second order systems 

utt ^ Po{D)u + Pi{D)ut (5) 

where 

s s 

j.k=i j=i 
We calculate simple wave solutions. Introducing 

gives us 

Utt = -\u;\^Poiij')u + i\u;\P,{u;')ut. (6) 

We have 

Lemma 1. A necessary condition for well posedness is that, for alluj' , the eigen- 
values k of 

[-k'^I + Pi(uj')k + Po{J)]a = Q (7) 

are real. 

Proof. If k{u>'),a{u') is a solution of (Q) then —k{ijj'),a{u}') is a solution if we 
replace uj' by —lo'. Since the solutions of (||) are of the form e'l"^''^" ^*a(a;') we 
only avoid catastrophic growth if k is real. 

If Pi = then becomes 

[-k^I + Pf){uj')]a = Q 

and Lemma 1 reduces to 

Lemma 2. If Pi ^ Q then a necessary condition for well posedness is that the 
eigenvalues of Pq{lu') are real and nonnegative. 

We can write as a first order system by introducing a new variable v with 
Ut = i\uj\v. We obtain 

(l).-M(-'r'T')(a)-"-lK«) 

The eigenvalues of P are determined by (|^). Using this reduction, we can define 
what we mean by strongly and weakly hyperbolic (second order) systems. 

Definition 3. The system is called strongly hyperbolic if for all uj' the eigen- 
values of P are real and there is a uniformly linearly independent (in lu') complete 
set of eigenvectors. 
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For strongly hyperbolic systems one can again develop a rather complete 
theory for local existence of quasi-linear systems. In particular lower order terms 

s 

Qu = BjDjU + B^ut + Cm 

do not destroy the well posedness of the problem. 
If Pi = we have 

Theorem 4. Assume that Pi = 0. The system is strongly hyperbolic if and only 
if the eigenvalues of Po{uj') are strictly positive and Pq{lu') has a complete set of 
eigenvectors which is uniformly ( in uj' ) linearly independent. 

Proof. Notice that when Pi = 0, any eigenvector of P with eigenvalue A_, (a;') is 
of the form 

(9) 



where the splitting corresponds to the block structure of P. Moreover, for each 
eigenvector of the form (||) there is another eigenvector 



with eigenvalue —Xj(uj'), which is linearly independent from the first if and only 
if Xj 0. 

Now, it is easy to check that a set of eigenvectors 

with real eigenvalues {Xj, — A^}, j = 1, 2, . . . , ti is a set of 2n uniformly linearly 
independent (in w') eigenvectors if and only if the set {aj{uj')}, j = 1,2,. ..n, 
is a set of uniformly linearly independent (in lo') eigenvectors of Pq[uj') with 
positive eigenvalues X'j{uj') > 0. This proves the theorem. 

Definition 4. The system (^ is called weakly hyperbolic if for all iv' the eigen- 
values of P are real. 

In particular we have 

Lemma 3. If Pi = and Pq{lu') has zero as an eigenvalue then the system is 
not strongly hyperbolic. It is weakly hyperbolic if the eigenvalues are real and non 
negative. 

As in the case of first order systems, weakly hyperbolic systems can have 
catastrophic exponential growth when adding lower order terms or considering 
variable coefficients. As example we consider 

Utt = aUxX + Uyy + bUx + CUy 
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The problem is strongly hyperbolic if a > and weakly hyperbolic if a = and 
there is catastrophic exponential growth if a = and & ^ 0. 

In the next section we will show that the ADM equations, linearized around 
flat, are only weakly hyperbolic for constant lapse function and shift vector. We 
can transform it into a strongly hyperbolic system if we choose the lapse function 
proportional to the trace of the solution. However, such a choice might introduce 
singularities. 

Consider, for example, 

uu = axUx + a^x- 

If a = a{x, t) is a given hmction, then the equation is weakly hyperbolic. If we 
choose a = u, we obtain 

Utt = [Uxf + Uxx- 

Now the equation is strongly hyperbolic but we will, in general, encounter sin- 
gularities due to the lower order term. 

3 Second Order Initial Value Formulations For General 
Relativity 

Our starting point are the ADM equations for General Relativity. The 3- 
metric induced on the spacelike 3-surfaces t — constant is denoted by all 
latin indices run over 1,2,3. 

From start we fix the shift vector equal to zero but keep a general lapse 
function a. Using the ADM equation for to eliminate the extrinsic curvature 
from the other ADM equation we get a second order evolution equation for 7^^- 

dt'-iij = a^i^"^{didmli] + didjjini - dida„ij - djOa^i^j + 2adidja 

+ lot, (10) 

where all derivatives are partial derivatives with respect to time and the coor- 
dinates on the t = constant 3-surfaces. Here and below "lot" stands for "lower 
order terms", that is functions of a, jij and their first derivatives. For the pur- 
pose of this paper, it is enough to say that all lower order terms are quadratic 
in first order derivatives. 

We have to consider the constraint equations. The momentun constraint is 

7'''at(9,7,/-9,7,i) +lot = 0, 

while the Hamiltonian constraint is 

7*^7'™(a,aajm - d.djjiJj + lot = 0. 

We now linearize around a flat solution (Minkowski spacetime) in Cartesian 
coordinates, that is we make 

7y =%+e/i,j- and 7'™ = 5'" - £ -I- ©(e^) (n) 
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and keep terms linear in e. Our new variable is hij. The constraint equations 
become 

dt{d,H~5''"'dih,„}j =0 

and 

5'^5'■'^^,^lh,ra-^H = Q 

where H = tr(/iij) = 5"^^ hij and A is the Laplacian in Ti? . Both linearized 
constraint equations are satisfied if 

a,iI-5'"9,/i„„ = 0. (12) 

Before linearizing (^0|) we have to make a choice of lapse. On the one hand 
the simplest possible choice a = 1 gives, after linearization, 



so that (12) implies 

d^hrj = Ah^j - ^^^JH. (13) 



On the other hand, choosing 

k k 
a = gtr(7y) = -(5'"7/™, 

gives after linearization and using ( p^ ) 



d^h,j = e (Ah,j - iftaji?) . (14) 

We define u — (ft-n, /i22, ^33, /ii2, /ii3, /i23)' to analyze the hyperbolicity of 
equations (|l^) and (p^. 

Thus the matrix Po{uj'), as defined in Section 2.2, of the system asociated 
to ( p^ ) has constant eigenvalues: fc^ with multiplicity five, and (2/3)fc^ with 
multipficity one. Also, this matrix has a uniformly linearly independent complete 
system of eigenvectors. Then according to Theorem 4, the system is strongly 
hyperbolic for any k ^ 0. < 1 gives a system with characteristic speeds 
smaller or equal than one, while |fc| > 1 would be an "unphysical" (though 
strongly hyperbolic) system with characteristic speeds higher than one. 

The matrix Po{i^') of the system asociated to (^3|) is uniformly diagonable 
but its eigenvalues are: with multiplicity one and 1 with multiplicity five. Thus, 
according to Lemma 3, equation (|l3| ) is only weakly hyperbolic. 

Another possibility to "cure" equation (|l^) is the usual choice of coupling 
the lapse function to the determinant of the 3-metric instead of coupling it to 
the trace as we have done above. This also leads to a strongly hyperbolic system. 
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4 Difference Approximations 

Consider the simple model problem 

Ut — Ux, —oo < X < oo, t > (15) 

with initial data 

uix,0)^f{x). (16) 

We are interested in solutions which are 27r-periodic in space. We want to solve 
the above problem by difference approximation. For that reason we introduce a 
gridlength h — 2Tr/N, N a natural number, gridpoints Xi, = vh and gridfunctions 
Uv{t) = u{xu, t). We also introduce the usual centered difference operators by 



jdx? ~ D+D^Ui, — {u^+i — 2uy + Uy-i)/h^ . 
Then we approximate (|l5|),([l^) by 

{u,)t = Dou,, j. = 0,l,2,...,iV-l (17) 
with periodic boundary conditions 

Ui^{t) ^ Ui,+N[t) (18) 

and initial conditions 

«.(0) = (19) 

(0) ^ (111) represents a system of ordinary differential equations which we solve 
with help of a standard ODE solver like the usual Runge-Kutta method. 
We want to discuss the accuracy of the approximation. We assume that 

M 

f(x)= e^"'/M' M<NI2. 

Then we can expand both the solutions of (p^,(p^) and (p7|)-(pg|) into Fourier 
polynomials 

M M 

u{x,t)= e"^''u{uj,t), u^{t)= Y e"^'""Ht^,t) (20) 

uj=~M uj = -M 

with 

u{iu,0) = Hiu,0)^ fiu). 
We introduce (|20| ) into (|l^) and respectively. Since 

oe /ox = iLoe and Df)e — e , 
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we obtain, for every frequency, 

Ut(UJ,t) = iLL!u[iLj,t), Ut[LJ,t) = iauju[uj,t), a— — . 

OJIl 

Therefore, 

uiLO, t) = e'^'fiu;), iL{cu, t) = e^""*/(^)- 
Thus there is a phase error 

e = (1 — a)Lot. 

Also, the wave speed for the difference approximation depends on the frequency, 
i.e., there is dispersion for the difference approximation but not for the differ- 
ential equation. This causes lots of difficulties if the solution is not properly 
resolved, i.e., there are not enough points/wavelength. 

Instead of the above second order method, one can use higher order methods. 
This results in a remarkable improvement of the accuracy. In Table |l| we give 
the number of points/wavelength so that the numeric solution has a phase-error 
of 10% or 1% after calculating during q time periods with methods of different 
order. 



Table 1. Points/ wavelength 



e 


2nd Order 


4th Order 


6th Order 


10% 


20 






1% 


64 


13 gi/^ 





P. Huebner and J. Thornburg using fourth order accurate methods, 
have demonstrated the improved accuracy for the Einstein equations. 

If we calculate with N points, then the solution consists of M ^ N/2 sim- 
ple waves. Most of them have large phase errors. Therefore, the approxima- 
tion is only useful if the Fourier expansion of the analytic solution converges 
rapidly. In particular, the part of the solution consisting of those waves with few 
points/wavelength travels in the wrong direction. 

Consider 

with highly oscillatory initial data 

u,y(0) = {—lYgu^ g smooth function. 
Introduce a new variable by Ui, = {—\Ywu- Then w solves 
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w{x,0) = g{x). 

Thus u^{t) represents a highly oscillatory wave which travels in the "wrong" 
direction. The usual way to control Uv{t) is to add an artificial viscosity term. 
As an example, consider 

Ut = xux, — 1 < a; < 1, 

with boundary condition 

u{-l,t) = -l, = 

and the initial data 

7r 

u{x, 0) = — cos — (x + 1). 

The solution forms an internal layer at a; = where the gradient becomes larger 
and larger. If we use the approximation 

Uut = x^Df)U^, 

then there will be a highly oscillatory wave traveling out of the layer. 
We approximate the wave equation 

Utt = Uxx 

by 

Utt = D+D-u. 

For the same level of accuracy, we need only half the number of points/wave- 
length. Also, there are no spurious waves which travel in the wrong direction. 



5 Constraints 

Using an example from fluid dynamics we want to demonstrate some of the 
problems which can arise when solving equations with constraints. Consider the 
Stokes problem 

Ut+Px = Aw, (21) 

Vt+Py =Av, (22) 

d=:Ux+Vy = 0, (23) 

in some domain x (0. T) with bomidary conditions 



Un = for (x, y) € dQ, 



0<t<T. 
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Here u, v denote the velocity components, u„ the normal component and p the 
pressure. Differentiating the first equation with respect to x and the second with 
respect to y gives us, using the divergence relation d — 0, 

Ap = 0. (24) 

We solve (§l|),(|2|) and (|2|) and think of Q as the constraint. 
One could be tempted to use, for p, 

p — po, {x,y) G dQ, pq given function, (25) 

as boundary condition. However, then we would, in general, not preserve the 
constraint d = 0. 

Differentiating (|2|) with respect to x and (|2^ ) with respect to y and using 
( p4[ ) gives us an equation for the divergence d, 

dt = Ad. 

By assumption, o? = for t = but we cannot guarantee that d = at later 
times if we use the boundary conditions (p5|). We must use 

d^O for {x,y) edQ, t > 0, 

as boundary condition and we cannot give p. 

Let Ah = Dj^xD-x + D^yD-y. A typical difference approximation is given 

by 

ut + DqxP = Ahu, (26) 
vt + Doyp = Ahi, (27) 
AhP = 0. (28) 

For the discrete divergence dh = Dq^u + D^yV we then obtain 

dht =-{Dlx + Dly)p + Ahdh, (29) 
AhP = 0. (30) 

The difficulty here is that 

A, ^ Dlx + Dly 
and therefore divergence is created. Instead of ( p8|) one can use 

AhP ~ ctdh, a » 1 constant. 

Then we can write ( p9| ) as 

dht + [DIx + Dly - Ah)p + adh = Ahdh- 

The damping term adh keeps dh under control. 
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